#include "includesT.h"
#include "pocsT.h"


using namespace std;


struct thread_data
{
	int J;
	int K;
	int N;
	ScanData *S;
};

struct thread_data thread_data_array[NUM_THREADS];

void *IterateLine(void *threadarg){
	struct thread_data *data;
	data = (struct thread_data *) threadarg;
	int j = data->J;
	int k = data->K;
	int n = data->N;
	ScanData *s= data->S;
	int i;
	for(i=1;i<=n;i++){
		P1hilbert(*s,j,k);
		P2compact(*s,j,k);
		P3apriori(*s,j,k);
		P4fproject(*s,j,k);
		P5nonneg(*s,j,k);
	}

	pthread_exit((void*) j);
}

void ThreadPocs(ScanData & s){

	int numit = 55;
	int rc;
	int j,k;
	/* Run each line in horizontal planes in parallel, one after the other*/
	//for(k=0;k<s.Nz;k++){
	for(k=255;k<=255;k++){
		/* Initialize and set thread detached attribute */
		void *status;
		pthread_t thread[NUM_THREADS];
		pthread_attr_t attr;
		pthread_attr_init(&attr);
		pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
		cout << "Iterating line" << k << endl;
		for (j=0;j<s.Ny;j++){
			if(s.Lmask[j][k]){
				thread_data_array[j].S = &s;
				thread_data_array[j].N = numit;
				thread_data_array[j].J = j;
				thread_data_array[j].K = k;
				rc = pthread_create(&thread[j], &attr, IterateLine, (void *) &thread_data_array[j]);
				if (rc) {
					cout<< "ERROR; return code from pthread_create() is "<< rc <<endl;
					exit(-1);
				}
			}
		}
		pthread_attr_destroy(&attr);
		for(j=0;j<s.Ny;j++) {
			if(s.Lmask[j][k]){
				rc = pthread_join(thread[j], &status);
				if (rc) {
					cout<< "ERROR; return code from pthread_join() is "<< rc <<endl;
					exit(-1);
				}
			}
		}
	}
}

int main (int argc, char* argv[]) 
{
	time_t t0, t1,t2;
	int i,j,k,p,n,q,P;
	bool initalguess = false;
	Float x=0,y=0,z=0;
	double dif;
	// Nu, umin, umax, Nv, vmin, vmax, Nt, tmin, tmax, R, D, Nx, Ny, Nz, zmax, xmax 
	cout << "Constructing ScanData"<< endl;
	ScanData s = ScanData(512, -255.5, 255.5, 512, -255.5, 255.5, 720, 3*PI/2, 3*PI/2, 416, 832, 512, 512, 512, 100, 100); // Non-truncated
	//ScanData s = ScanData(116, -116, 116, 256, -256, 256, 360, 3*PI/2, 3*PI/2, 416, 832, 512, 512, 2, 1, 128); // Truncated
	cout << "Reading data" <<endl;

	// Compact support data test
	Float a=69*1.05,b=91*1.05, c=91*1.05;
	for (i=0; i<s.Nx; i++){
		x = s.xmin + i*s.dx;
		for (j=0;j<s.Ny; j++){
			y = s.ymin + j*s.dy;
			for(k=0;k<s.Nz;++k){
				if (x*x/a/a + y*y/b/b + z*z/c/c < 1.) s.CSmask[i][j][k] = true;
			}
		}
	}
	
	
	ifstream file("../../../data/cbct/test2.prj", ios::binary | ios::ate | ios::in); // Non-truncated
	//ifstream file("../../../data/cbct/Data_trunc.bin", ios::binary | ios::ate | ios::in); // Truncated
	if (file.is_open()) {
		long size = file.tellg();
		file.seekg(0, ios::beg); 
		float* unsorted;
		unsorted = new float[size];

  		file.read((char*)unsorted, size);
  		file.close();
		unsigned long N =0;
		for (p=0;p<s.Nt;++p){
			if(p<s.Nt/4){ P = p + 3*s.Nt/4;
			}else{ P = p - s.Nt/4;}
			for(q=0; q<s.Nv;++q){
				for(n=0; n<s.Nu;++n){
					s.data[P][n][q] = unsorted[N];
					N++;
				}
			}
		}
		delete[] unsorted;
	}else{cout<< "Unable to open data file";}

	
	file.open("../../../data/cbct/APmask_512_circle.bin",ios::binary | ios::ate | ios::in);
	if (file.is_open()) {
		long size = file.tellg();
		file.seekg(0, ios::beg); 
		int* unsorted;
		unsorted = new int[size];

  		file.read((char*)unsorted, size);
  		file.close();
		unsigned long N =0;

		for (k=0;k<s.Nz;++k){
			for(j=0; j<s.Ny;++j){
				for(i=0; i<s.Nx;++i){
					s.APmask[i][j][k] = (bool) unsorted[N];
					N++;
				}
			}
		}
		delete[] unsorted;
	}else{cout<< "Unable to open APmask file";}

	file.open("../../../data/cbct/APdata_512.bin",ios::binary | ios::ate | ios::in);
	if (file.is_open()) {
		long size = file.tellg();
		file.seekg(0, ios::beg); 
		float* unsorted;
		unsorted = new float[size];

  		file.read((char*)unsorted, size);
  		file.close();
		unsigned long N =0;
		for (k=0;k<s.Nz;++k){
			for(j=0; j<s.Ny;++j){
				for(i=0; i<s.Nx;++i){
					s.APdata[i][j][k] = unsorted[N];
					N++;
				}
			}
		}
		delete[] unsorted;
	}else{cout<< "Unable to open APdata file";}

	if(initalguess){
		file.open("recon.bin",ios::binary | ios::ate | ios::in);
		if (file.is_open()) {
			long size = file.tellg();
			file.seekg(0, ios::beg); 
			float* unsorted;
			unsorted = new float[size];

			file.read((char*)unsorted, size);
			file.close();
			unsigned long N =0;
			for (k=0;k<s.Nz;++k){
				for(j=0; j<s.Ny;++j){
					for(i=0; i<s.Nx;++i){
						s.recon[i][j][k] = unsorted[N];
						N++;
					}
				}
			}
			delete[] unsorted;
		}else{cout<< "Unable to open Initial guess file";}
	}

	time(&t0);
	cout << "Backprojecting data" << endl;
	prepare(s);
	time(&t1);
	cout << "Preparation time: "<< difftime(t1,t0)<<" sec"<<endl;
	cout << "Iterating" << endl;
	ThreadPocs(s);
	time(&t2);
	cout << "Time to iterate: "<< difftime(t2,t1)<<" sec"<<endl;
	cout<< "Writing to file" << endl;

	ofstream outfile("recon.bin", ios::binary | ios::trunc);
	if (outfile.is_open()) {
		Float* unsorted;
		long size = sizeof(Float)*s.Nx*s.Ny*s.Nz;
		unsorted = new Float[size];
		unsigned long N=0;
		
		for (k=0;k<s.Nz;++k){
			for(j=0; j<s.Nx;++j){
				for(i=0; i<s.Ny;++i){
					unsorted[N] = s.recon[i][j][k];
					N++;
				}
			}
		}
  		outfile.write((char*)unsorted, size);
  		outfile.close();

		delete[] unsorted;
	}else{cout<< "Unable to open data file";}

	
	return 0;
}
